function [robot_his,cont_his] = runScenarioOnce(r_num,duration)
    if nargin < 2
        duration = 500;
    end
% given the number of robots, simulate its serving process for one fixed hotspot center
    warning('off'); close all; 
    % filed limits [x_min x_max y_min y_max z_min z_max]
    field_limits = [0 100 0 100 0 20];
    robot = [ rand(r_num,1)*field_limits(2), rand(r_num,1)*field_limits(4), rand(r_num,1)*field_limits(6) ];    
    % plot of trajectory
    
    % discretizing the x-y plane, 
    discretization_num = 50; mesh_size = field_limits(2)/discretization_num;
    x_space = linspace(field_limits(1),field_limits(2),discretization_num);
    y_space = linspace(field_limits(3),field_limits(4),discretization_num);
    [X_points,Y_points] = meshgrid(x_space,y_space);
    points_xy = [reshape(X_points,discretization_num^2,1) reshape(Y_points,discretization_num^2,1)];

    % density calculation
    centers = [50 50]; c_index = 1;
    d_center = centers(c_index,:);
    R = densityFunc(d_center,points_xy);
    sensor_para = 20;

    % variables preparation 
    robot_his = cell(duration,1);     robot_his{1} = robot; cont_his = cell(duration,1); 
    robot1 = robot;
    for k = 1:duration
        control = zeros(size(robot1)); 
        for i = 1:r_num
            % refine the robot set of multual communication
            s_i = robot(i,:); 
            rel_dis = vecnorm(s_i(1:2) - robot(:,1:2),2,2); 
            r = coverageModel(robot(:,end),sensor_para,'reciprocal');
            robot_i = robot; r1 = rel_dis < r + r(i); r1(rel_dis == 0) = 1;
            robot_i = robot_i .* r1;
            [r,dr_dh] = coverageModel(robot_i(:,end),sensor_para,'reciprocal'); r = r.*r1; dr_dh = dr_dh.*r1;
            % compute the gradient
            delta_u_i = computeGradient(i,robot_i,r,dr_dh,points_xy,mesh_size,R,d_center);
            direction = delta_u_i/norm(delta_u_i,2);
            control(i,:) = 0.1* direction';
        end
        robot1 = robot + control;
        k/duration;
        robot_his{k+1} = robot1;        robot = robot1; robot(:,3) = abs(robot(:,3));
        cont_his{k} = control;
    end    
end

